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We study numerically statistical distributions of sums of orbit coordinates, viewed as indepen- 
, ^ dent random variables in the spirit of the Central Limit Theorem, in weakly chaotic regimes 

rN ■ associated with the excitation of the first (fc = 1) and last {k = N) linear normal modes of 

I the Fermi-Pasta-Ulam-a system under fixed boundary conditions. We show that at low energies 

{E = 0.19), when the fc = 1 linear mode is excited, chaotic diffusion occurs characterized by dis- 
tributions that are well approximated for long times (t > 10^) by a g-Gaussian Quasi-Stationary 
State (QSS) with q « 1.4. On the other hand, when the k — N mode is excited at the same 
energy, diffusive phenomena are absent and the motion is quasi-periodic. In fact, as the energy 
increases to _E = 0.3, the distributions in the former case pass through shorter g-Gaussian states 
and tend rapidly to a Gaussian (i.e. q — )• 1) where equipartition sets in, while in the latter we 
need to reach to i? = 4 to see a sudden transition to Gaussian statistics, without any passage 
through an intermediate QSS. This may be explained by different energy localization proper- 
ties and recurrence phenomena in the two cases, supporting the view that when the energy is 
placed in the first mode weak chaos and "sticky" dynamics lead to a more gradual process of 
energy sharing, while strong chaos and equipartition appear abruptly when only the last mode 
is initially excited. 

Keywords: Non-Extensive Statistical Mechanics, g-Gaussian distributions, "edge of chaos", weak 
chaos, quasi-stationary states, multi-dimensional Hamiltonian systems 
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1. Introduction 

Probability density functions (pdfs) of chaotic trajectories of dynamical systems have been studied for 



Anosovl. 


19671: lArnold k Avez 


198C. 198 


2I: Eckmann & Ruelle 



19851 : iKatok et al 



Pesinl . Il976l . Il977l : iRuelld . Il979 



198,5l ]. The fundamental questi 



ion 



in this regard concerns the existence of an appropriate invariant probability density (or ergodic measure) 
which characterizes chaotic motion in phase space regions where solutions generically exhibit exponential 
divergence of nearby trajectories. If such an invariant measure can be found for almost all initial conditions 
(i.e. except for a set of measure zero), one has a firm basis for studying the system from a Statistical 
Mechanics point of view. 

If this invariant measure turns out to be a continuous and sufficiently smooth function of the phase 
space coordinates, one can invoke the Boltzmann-Gibbs (BG) microcanonical ensemble and attempt to 
evaluate all relevant quantities of equilibrium Statistical Mechanics, like partition function, free energy, 
entropy, etc. On the other hand, if the measure is absolutely continuous (as e.g. in the case of the so- 
called Axiom A dynamical systems), one might still be able to use the formalism of modern ergodic 
theory and SRB measures (after Sinai, Ruelle and Bowen) to study the statistical properties of the model 
Eckmann &: Ruelld . [l985l | . 



Thus, viewing the values of one (or a linear combination) of components of a chaotic solution at discrete 
times tn, = 1, . . . , A/" as realizations of M independent and identically distributed random variable s X„ , 
and calculating the distribution of their sums, in the sense of the Central Limit Theorem (CLT) [ Ricd . 
1995l |. one expects to find a Gaussian pdf, whose mean and variance are those of the A„'s. This is indeed 
what happens for many chaotic dynamical systems studied to date which are ergodic, i.e. almost all their 
orbits (except for a set of measure zero) pass arbitrarily close to any point of the constant energy manifold, 
after sufficiently long times. In these cases, at least one Lyapunov exponent is positive, stable periodic 
orbits are absent and the constant energy manifold is uniformly covered by chaotic orbits, for all but a 
(Lebesgue) measure zero set of initial conditions. 

It is also important, however, to study "small size" chaotic regions of Hamiltonian systems, at energies 
where the maximum (positive) Lyapunov exponent is "small" and stable periodic orbits exist with islands 
of invariant tori and sets of cantori, which occupy a positive measure subset of the energy manifold. In such 
regimes of "weak chaos" , a great number of orbits stick for long times to the boundaries of these islands and 
chaotic trajectories diffuse slowly through multiple conne cted regions in a highly non-uniform way [Aizaw 
1984 : IChirikov &: Shepelyanskyl . 11984 : iMeiss &: Ottl.ll986 1. Such examples occur in many physically realistic 



syste r ns studied in the cu rrent literature (see e.g. [Skokos et al. 
20091 : ISkokos et ajl l2009l ]l 



20081 : iFlach et a.ll l2009l : I.Tohansson 



In this paper, we plan to focus on such "weakly chaotic" regimes of the first and last normal modes 
of the Fermi-Pasta-Ulam (FPU)-a 1-dimensional particle chain under fixed boundary conditions. We thus 
demonstrate, by means of numerical experiments according to the CLT, that pdfs of sums of orbit compo- 
nents at low energies do not rapidly converg e to a Gaussi an, but are well approximated for long integration 
times by the so-called g-Gaussian function Tsallis . [2009l | 



P{s) = dexp (-/3s^) = d 



1 - (1 - q)l3s'' 



1 

l-q 



(1) 



where the index satisfies 1 < g < 3, /3 is an arbitrary parameter and d is a normalization constant. At 
longer times, of course, weakly chaotic orbits are expected to eventually escape to larger chaotic seas, 
where obstruction by islands and cantori is less dominant and the dynamics is more uniformly ergodic. 
This transition is associated with energy equipartition among all linear modes and is characterized by 
g-Gaussian pdfs ([T]) whose index decreases to q = 1, representing the limit at which the pdf becomes a 
Gaussian. 

Thus, in our models, g-Gaussian distributions represent quasi- stationary states (QSS) that are often 
very long-lived, especially inside "thin" chaotic layers, where the Lyapunov exponents are small. This 
suggests that it might be useful to study these pdfs (and their associated q values) for sufficiently long 
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times to derive information about the long term dynamics of the QSS and their connection to energy 
equipartition among all degrees of freedom. 

In what follows, after introducing our methodology in Sec. [21 we begin in Sec. [3] by describing the 
system under study and some of its main properties. Next, in Sec. [Hwe study two special solutions of the 
Fermi-Pasta-Ulam-a (FPU-a) 1-dimensional chain of = 31 particles under fixed boundary conditions. 
These solutio ns correspond to: (a ) The excitation of the k = 1 linear mode, which has been related to the 
FPU paradox Fermi et al\ . \l95d! \ and (b) the excitation of the last {k = 31) linear mode, at energies where 
the nonlinear continuations (a > 0) of these modes have turned unstable. In case (a), we find that chaotic 
diffusion phenomena (related to the breakdown of the so-called FPU recurrences) are characterized by 
g-Gaussian pdfs with q decreasing to 1 as energy increases and equipartition sets in at earlier and earlier 
times. On the other hand, case (b) shows persistent quasi-periodic behavior until, it suddenly becomes 
strongly chaotic and equipartition abruptly occurs. These results may be explained by examining the 
exponential localization properties of the averaged modal energies of the two solutions in A;-space. 



2. Statistical distributions of chaotic QSS and their computation 

The problem we investigate is described by an autonomous N degree of freedom Hamiltonian function of 
the form 



H = H{q{t),p{t)) = H{qi{t), . . . , <7^(t),pi(t), . . . ,pN{t)) = E 



(2) 



where {qk{t),Pk{t)) are respectively the positions and momenta representing the solution in phase space at 
time t. As is well-known, these solutions can be periodic, quasi-periodic or chaotic depending on the initial 
condition (g(0),p(0)) and the value of the total energy E. What we wish to study h ere is the statistics of su ch 



systems in regimes of "w e akly" chaotic motion, where the Lyapunov exponents Benettin et al . 1980al jbl: 



Eckmann fc Buellel . Il98,4 ISkokod . l2ninl | are positive but very small. Such situations often arise when one 
considers solutions which diffuse slowly in thin chaotic layers and wander through a complicated network 
of higher order resonances, often stic king for very l ong times to the boundaries of islands const i tuting 
the so-called "edge of chaos" regime [Aizawal . Il984l : IChirikov fc Shepelyanskyl . Il984l : iMeiss fc Ottl . Il986l : 
Tsalligl . I2nn9l |. 

At this point one can ask the following questions: How long do these "weakly" chaotic states last? 
Assuming they are quasi-stationary, can we describe them statistically by following some of their chaotic 
orbits? What type of distributions characterize these QSS and how could one connect them to the actual 
dynamics of the solutions of the corresponding Hamiltonian system? 

Following an approach inspired by the well-known Central Limit Theorem Ricd . 1995l |. we shall use 
the solutions of Hamilton's equations of motion 



dt 



dH 



dPk 
dt 



dH 

dqk 



A; = 1, 



(3) 



to construct distributions of suitably rescaled sums of M values of a generic observable rji = r]{ti) {i = 
1, . . . , M) which depends linearly on the components of the solution. If these are viewed as independent 
and identically distributed random variables (in the limit M — t- oo), we may evaluate their sum 



5' 
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M 



M 



1=1 



(4) 



for j = 1, . . . , Nic different initial conditions. Thus, we can study the statistical behavior of the variables 



S^j , centered about their mean value {S^^j) = 77- Sj=i "^iLi vl''^ and rescaled by their standard deviation 



Jj) 



S3) 
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M N,, M ^ 
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where 



2 



(6) 



Plotting the normaHzed histogram of the probabihties P{s^J) as a function of s^j , we then compare our 
pdfs with a g'-Gaussian function of the form 



P{sf) = dexp^(-/34)2) ^ d 



1 



(i)2 



1 

1-9 



(7) 



where q is the so-called entropic index, /3 is a free parameter and d a normalization constant [Xsallid . 
2009i]. Note that, in the limit g — )■ 1, ([7]) tends to the well-known Gaussian distribution, i.e. expg(— — >■ 
exp(— 

It can be shown that the g-Gaussian distribution ([7|) is normalized when 



/3 = d^/^- 



3-9 
2(9-1) 



(9-i)^r(^ 



(8) 



where T is the Euler T function j Tsallis . 2009f |. Clearly, ^ shows that th e allowe d values of q are 1 < q < 3. 
The index q appearing in ([7]) is connected with the Tsallis entropy Tsallis . [2009l] 



Sn k- 



yw q 

a-1 



w 



with ^^Pi = 1 



(9) 



where i = 1, . . . ,W counts the microstates of the system, each occurring with a probability pi and k is the 
so-called Boltzmann universal constant. Just as the Gaussian distribution represents an extremal of the 
BG entropy S-qq = Si = k Yl^i Pi liiPii so is the (/-Gaussian ([7]) derived by optimizing the Tsallis entropy 
dH]) under appropriate constraints. Systems characterized by the Tsallis entropy are said to lie at the "edge 
of chaos" and are significant lv different than BG systems, in the se nse that their entropy is non-additive 
and generally non- extensive |Tsallisl . I2nn9l : iTsallis fc Tirnaklil . l2ninl | . 

We now describe the numerical approach we follow to calculate the above pdfs. First, we specify an 
observable denoted by r]{t) as one (or a linear combination) of the components of the position vector q{t) of 
a chaotic solution of located initially at (g(0),p(0)). Assuming that the orbit visits all parts of a QSS 
during the integration interval < t < tf, we divide tj into Nic equally spaced, consecutive time windows, 
which are long enough to contain a significant part of the orbit. Next, we subdivide each such window into 
a number M of equally spaced subintervals and calculate the sum S^^J of the values of the observable ri{t) 
at the right edges of these subintervals (see eq. ([1])). 

In this way, we treat the beginning of every time window as a new initial condition and repeat this 
process Nic times, to obtain as many sums as required for reliable statistics. Consequently, at the end of 
the integration, we compute the average and standard deviation of the sums dU, evaluate the Nic rescaled 
quantities and plot the histogram P{s^^) of their distribution. 

As we shall see in the next sections, in regions of "weak chaos" these distributions are well- fitted by 
a g-Gaussian of the form ^ for fairly long time intervals. However, for longer times (or higher energies), 
as the orbit diffuses in domains of "strong chaos", a Gaussian pdf (with g ~ 1) is expected to describe the 
dynamics. 



3. The Fermi-Pasta-Ulam-a model 

The well-known FPU model is a Hamiltonian lattice of N particles (degrees of freedom) with a potential 
that depends only on nearest neighbor particle interactions. Its name was given aft er the pioneering work 
by E. Fermi, J. Pasta, S. Ulam and M. Tsingou in the 1950's [Fermi et all . \l95^ . which demonstrated 
the lack of thermalization of the system in contrast to what is expected by BG Statistical Mechanics. In 
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particular, they reported the absence of energy transfer to all normal modes of the system even after very 
long integration times, due to the so-called FPU recurrences, where the energy returns from neighboring 
normal modes to the first one that was initially excited. This phenomenon, termed "the FPU paradox", 
still remains an open problem, although much progress has been made to date towards understanding its 
dynamical behavior in the thermodynamic limit of E and N becoming arbitrarily large, with E/N =const. 

More specifically, in the present work we deal with the FPU-a model of identical particles on a 
1-dimensional chain described by the Hamiltonian 

n=l n=0 n=0 

where Qn is the displacement of the nth particle from its equilibrium position, pn is the corresponding 
momentum, a > is a real constant and E is the fixed energy of the system. Moreo ver, we shall consider 
here only fixed boundary conditions, i.e. qo = qn+i = following [Fermi et all Il955 |. 
As is already known, under the canonical normal mode transformation 

N 



9n=A/^^I]QfcSin(^A;n^^j , (11) 



Pn 



N 



2 . /, vr 



-VPfcsin A;n ,n = l,...,iV, (12) 

k=l 



system (|10p can be studied following the time evolution of the harmonic normal modes of the linear case 
(i.e. Q = 0), whose energies Ek are given by 

E, = Pi + iolQl (13) 

where 

Wfc = 2sin ( — — -],k = l,...,N (14) 



^2(iV + l 

are the harmonic frequencies. Taking the time averaged harmonic spectra 



1 



Ek{t) = - / Ek{s)ds (15) 
t Jo 

one expects to find that the system reaches equipartition of its total energy among all its normal modes if 
Ek H> E/N Vfc as t ^ oo. 



4. Single site excitations in the FPU-q model 

We now proceed to study two different cases of single site excitations in the FPU-a system (see eq. pOj) ) 
to investigate how energy spreads between all normal modes. The first is the well-known case of the 
excitation of the first normal mode k = 1, which is the one studied originally by Fermi, Pasta and Ulam in 
Fermi et al. 1, [^55] and the second one concerns the excitation of the last mode k = N, which correspond 



to the longest and shortest wavelength of the spectrum respectively. As we discuss below, they turn out to 
be strikingly different in the way they lead to the excitation of the full spectrum of E^, k = 1, . . . ,N. 

For the numerical integration of our trajector ies, we typically use throughout the paper Yoshida's 
fourth order symplectic integrator Yoshidal . Il990l | with a fixed time step varying between 0.01 — 0.05 
depending on the particular energy E of Hamiltonian (jlOh . which always results in a relative energy error 
of the order of 10~^ — 10~^. For the computation of the Lyapunov exponents and the Generalized Alignment 



Index (GALI) [Skokos et al\ . l2007l | using the above symplectic integrator (or any sym plectic integrator in 
gener al) , we apply the methodology of the Tangent Dynamics Hamiltonian proposed in (Skokos &: Gerlachl . 
2010'] which is suitable for the evolution of deviation vectors in the tangent space of the orbit under study. 
Having thu s access to the devi ation vectors, we compute the Lyapunov exponents following the standard 
method of Benettin et a/. . [l980a >b.] and GALIj, i>2 and use a numerically faster implementation of the 
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GALIs whic h is suitable especially for mult i-dimensional systems, called the "Linear Dependence Index" 
(LDI) as in lAntonopoulos &: Bountis . 2006]. 



Following [Ponno et al 



we s t udy t he case of A'" = 31 particles of the FPU-a system pOj) with 



a = 0.33. In particular, in [Ponno et all |2010| | a detailed study of energy diffusion is carried out, from the 
first normal mode k = 1 which is initially excited, to all remaining modes of the spectrum and an estimate 
of the time needed to reach equipartition is presented. Here, we are interested in focusing on regimes 
closer to the so-called "edge of chaos" , i.e. phase space regions where diffusion pro c esses are extremely 
weak due to the " s tickin ess" phenomenon [Aizawa . 19841 : Chirikov fc Shepelyansky . 1984 : Meiss fc Ott . 
19861 : Skokos et al . 2008]. Aided by the GALI meth od and t he co mputation of Lyapunov exponents, we 
use concepts of Non-Extensive Statistical Mechanics Tsallisl . 120091 ] to show that g-Gaussian distributions 
can describe quite satisfactorily weakly chaotic dynamics in excitations of the k = 1 mode, while they 
identify a distinctly different sharp transition to strong chaos and equipartition in excitations of the k = 3 1 
mode. Our results are in agreement with similar investigations carried out in Antonopoulos et al . 2010l |. 
concerning QSS in multi-dimensional FPU-/3 models. 



4.1. Excitations of the first normal mode 

Let us start our study exciting the first normal mode only, i.e. /c = 1 in eqs. pT]) and (fT2]) . In Fig. [1^) we 
present in logarithmic scale the averaged energy spectra of the system at time t = 10^ for the energy E = 
0.19. As has already been pointed out in Ponno &: Bambusil . l2005l : Bambusi &: Ponnol . [200i : iPonno et all 
2010l ]. the system is in a m etastable state, where f ew (low k) modes share the total energy of the system. 



forming a "natural packet" Berchialla et all 120041 ] , while all the rest gain a small amount of it and remain 



for long times exponentially localized in normal mode space. Although, for these parameters and initial 
conditions, the first normal mode seems to dominate the moti on, a rather weak d iffusion takes place in the 
highest modes of the spectrum, as reported and quantified in [Ponno et aLl . 120101 ]. In Fig.[T)3) the evolution 
of the averaged harmonic energies E^/E of eq. ()13p with N = 31, a = 0.33 and E = 0.19 is presented in 
log - log scale. The lower modes are depicted by the darkest colors and the higher by the lighter ones. As 
is evident, a weak diffusion of energy from the lowest to the highest modes is revealed, resulting in a slow 
increase of the Ek values for large k. This suggests that we can expect energy equipartition after a long 
time interval ~ 10^ during which the maximal Lyapunov exponents decrease towards zero, saturating 
thereafter at values close to 10~^, as shown on panel c) of Fig. [TJ 

It is therefore important to study the behavior in the immediate neighborhood of this mode at E = 0.19, 
in terms of statistical distributions by the methods described in Sec. [2l In panel d) of Fig. [T] we plot the 
numerical (red curve), g-Gaussian (green curve) and Gaussian (blue curve) distributions that correspond 
to final integration time tf = 10^^ using Nic = 10^ time windows and M = 100 terms in the computation 
of the sums for ry = qi + q2- Our results show that this weakly diffusive motion towards equipartition is 
described by a pdf that is well-fitted by a g-Gaussian of the form ([7]), with q ~ 1.416, ^ 5 x 10^^ for 
tf = W^ and q ^ 1.342, ~ 2.41 x 10"^ for i/ = 10^^ (see Fig. [Hi)). 



4.2. Excitations of the last normal mode 

In system (jlOp of = 31 degrees of freedom, the last normal mode is that of = 31 and has the shortest 
wave length. Since its frequency u^i ~ 1.9975 ~ 2 lies outside the acoustic resonance regime, = 31 is 
more strongly localized in modal space than the k = 1 mode and thus it is interesting to study it in its 
own right. 

In Fig. [2^) we present the averaged energy spectra of the k = 31 mode at t = 10^ and in panel b) 
the time evolution of it s averaged energ i es EU E for E = 0.19 and a = 0.33. In Flach et al . 20051 . 2006 : 
Ivanchenko et aj\ . 20061 : Flach &: Ponno . 2008| the authors have predicted that this kind of trajectories, 
for values of the parameters as considered here, lie close to a q-breather and thus the energy spectra that 
appear in Fig. [2^) are approximately the same as those of the q-breather. The energy spectrum of the latter 
arises from the continuation of the simple periodic orbit that corresponds to mode fc = 31 of the linear 
system (see eq. (llOp for a = 0). Our results also demonstrate that the FPU trajectory arising from the 
excitation of the mode fc = 31 is close to a q-breather, since the hierarchy of the harmonic energies E^ of 
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Fig. 1. FPU-a system with A*' = 31, a = 0.33 and E = 0.19. The initial condition is such that all the energy E is given to 
the first normal mode = 1 (see eq. (|11|1 and (I12p ). In panel a) we plot the normalized averaged energy spectrum a.t tj = 10^, 
while panel b) shows the time evolution of the averaged energy spectra Ej^/E of eq. (|15[l . Panel c) presents the time evolution 
of the four biggest Lyapunov exponents. Panel d) is a plot (in lin - log scale) of the numerical data (red curve), g-Gaussian 
(green curve) and Gaussian (blue curve) distributions that correspond to final integration time = 10^^ using A^^^ = 10^ 
time windows and AI — 100 terms in the computation of the sums for the observable = gi + g2- In this case the numerical 
fitting with a g-Gaussian gives q « 1.342 with ~ 2.41 x 10~^. 



the q- breather, calculated by usi n g the proposition in IPenati fc Flachl . l2007l ] for a single mode excitation 



or m 



Christodoulidi et all \201Ui \ , [Christodoulidi et all l201ll | for multiple mode excitation, is almost the 



same as the tra.iector y shown in Fig. [2^). 

In particular, in 
that if we start with k 



Christodoulidi et all 1201 1[ | it is predicted, using the Poincare - Lindstedt method, 
31 at zeroth order, the modes ki = 2, k2 = N - 2 = 29, ks = 4, k4 = N - 4 = 27 
etc. will be successively excited at higher orders (for more details see the Appendix below). Knowing the 
sequence t hey follow, fci , fco , . . . , fc/y, we can identify these modes with those called "tail modes" in the 
literature [Penati k Flachl . boOTi : IPonno et all I2OI0I ] as having the lowest energy values. Thus, choosing 
these modes to be k 
dominant modes k 



1) • 



jv/2j • • • J one expects that the energy will diffuse from the strongly localized and 
, . . , kj\i/2-i to the tail modes. However, as we deduce from Fig. [2)3), there is no presence 
of energy diffusion to the tail modes of the system at least up to = 10®. 

This suggests that the localization phenomenon is very strong compared with that of the k = 1 case, 
on the same energy manifold E = 0.19. By computing the LDIj =GALIj, i = 2,..., 8 in Fig. [2h), we 
assert that the motion remains ordered at least up to t/ = 10® that we have checked numerically, in fully 
agreement with what we have found in panel b) of the same figure. The same finding is also justified by the 
fact that the four maximal Lyapunov exponents (not shown here) keep decreasing to zero up to the same 
final integration time. We conclude, therefore, that exciting only k = 31 yields a regular, quasi-periodic 
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t 



Fig. 2. FPU-a system with A'^ — 31, a — 0.33 and E = 0.19. The initial condition here is such that all the energy E is 
attributed to the last normal mode k = 31. Panel a) shows a plot of the normalized averaged energy spectrum at tj — 10^, 
exhibiting perfect exponential localization. In panel b) the time evolution of the averaged energy spectra Ej^/E of eq. (|15|l 
shows complete absence of diffusion, while in panel c) the time evolution of the LDI^ =GALIi, i = 2, . . . , 8 indicates that the 
motion is indeed quasi-periodic lying on a 2-dimensional torus. 



orbit, at least up to tj = 10^, apparently lying on a 2-dimensional torus, since GALI2 is the only index 
that remains nearly constant in this time interval. 

This striking difference in the time evolution of the harmonic energies in the above two cases is due to 
the acoustic resonance of k = 1 with a few nearby modes, i.e. uj^ — k ui. These modes are known to compose 
a "natural packet" and undergo an internal equipartition of the total energy. This natural packet has width 

1/4, 



equal to /i = a^/^ + l)i as predicted in Berchialla et all |2004| . l2005l : iBambusi Ponnol . I2OO6 



Ponno fc BambusJ - linn^l ]. For example, in Fig. [T]we obtain /x = 5.10279. 

By contrast, when we excite the fc = 31 normal mode, no diffusion through a metastable state is 
observed at energies E > 0.19, until about = 4 as we see in Fig. [3l In panels a), b) and c) of that figure, 
we show evidence of energy transfer from the dominant modes, i.e. A/^, 2, — 2, 4, — 4, 6, — 6, . . . to 
the tail modes with the lowest energy: 16, A^ — 16, ... , 30, A^ — 30. This transition occurs rather rapidly 
compared with the k = 1 excitations, since at times up to t/ = 10*^, energy equipartition between all 
modes has already taken place, as suggested by panel d) of Fig. [3l where we have plotted in lin - log scale 
the numerical (red curve) and Gaussian (blue curve) distributions for this integration time and observable 
T] = gi4 + qi5 + qiQ. These results clearly show that the sum distributions of the data are very close to 
Gaussians {q = 1), implying that by that time our orbit has entered a regime of strong chaos. 
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Fig. 3. FPU-Q system with = 31, a = 0.33 and E — 4. The initial condition here is such that all the energy is attributed 
to the last normal mode fc = 31. Panel a) plots the normalized averaged harmonic energy spectrum at tj = 10^ showing how 
the total energy has been distributed among the modes. Panel b) shows the time evolution of the averaged energy spectra 
Ef^/E of eq. (|15[l . while panel c) presents the time evolution of the four largest Lyapunov exponents. In panel d), plotting the 
numerical distribution (red curve) that corresponds to the final integration time = 10^, for A'^^ = 10^ time windows and 
M = 100 terms in the sums, we find for the observable rj — qi4 + gi5 + gig that the statistics is closely approximated by the 
Gaussian distribution {q = 1) denoted by the blue curve and hence the motion has entered a regime of strong chaos. 

5. Conclusions 

In this paper, we have numerically constructed pdfs of rescaled sums of observables derived from periodic 
orbits of the FPU-a Hamiltonian system. Assuming that these observables behave as independent random 
variables, we have attempted to determine their statistics in phase space regions of "weak chaos" in the 
context of the Central Limit Theorem. In particular, we focused our study on initial excitations of the 
first and last normal modes of the system and have examined the resulting dynamics from the point of 
view of energy diffusion processes leading to the formation of QSS, which eventually reach equipartition 
at sufficiently high energies and/or integration times. 

At energies where the periodic orbits have turned unstable, the motion is generally expected to evolve 
within weakly chaotic domains, where the associated pdfs are well approximated by (^-Gaussian distribu- 
tions {1 < q < 3) for long times. These pdfs frequently represent QSS and tend to Gaussians {q = 1) 
for longer integration times and/or higher energies, as the orbits diffuse into domains of "strong chaos" 
characterized by large (positive) Lyapunov exponents. 

In the case of the k = 1 mode excitation, the energy flow to higher modes, as predicted by the theory 
of FPU recurrences and the breakdown of exponential localization of the energies was found to be 
connected with the appearance of a QSS whose statistics is well-described by a g-Gaussian with q ^ 1.34, 
at a total energy oi E = 0.19. At higher energies, e.g. E = 0.3, this QSS is more short-lived and energy 
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equipartition occurs more rapidly, described by g-Gaussian pdfs which quickly converge to Gaussians with 
q = l. 

The excitation of the last normal mode, however, yields distinctly different dynamics: At low energy 
values, exponential localization of linear mode energies is very dominant so that no energy diffusion to 
higher modes is observed and no metastable states are formed, as the motion continues to oscillate quasi- 
periodically on low-dimensional tori. This persists until the total energy becomes w 4, where sum 
distributions rapidly converge to Gaussians and energy equipartition occurs abruptly as the motion enters 
suddenly a regime of strong chaos. 

Other authors have also studied similar QSS in conservative systems such as coupled sta ndard maps and 



the so -called HMF model from the point of view of Non-Extensive Statistical Mechanics [Baldovin et al. 



2nn4al lbl|. In these works, the growth of variance of a temperature-like quantity, expressed by the averaged 



"angular momentum" of the models, was calculated and interesting results were obtained, but not from 
the viewpoint of sum distributions. The connection of such metastable states to dynamical phenomena 
such as the appearanc e of chaotic breathers and the onset of equipartition has been analyzed in detail by 
Cretegny et al\ . 19981 ] in the vicinity of the 7r-mode in FPU-/3 systems with periodic boundary conditions. 



More recently, long-lived QSS approxim ated by q-Gaussian pdfs hav e been investigated in various multi- 
dimensional FPU-/? Hamiltonian systems in [Antonopoulos et al. . 20101 ]. under periodic and fixed boundary 
conditions. It was found that near unstable periodic orbits representing continuations of harmonic modes of 
the linear (/3 = 0) problem there exist weakly chaotic domains where pdfs of the g-Gaussian type describe 
QSS for surprising long time intervals. Beyond these intervals, g — )• 1 and the pdfs quickly converge to 
Gaussians, as the orbits enter regimes o f strong chaos and energ y equipartition occurs. Finally, certain very 
recent results on area-preserving maps [Ruiz-Lopez et al. . 20 IQ] are also in good agreement with what we 
have found here for multi-dimensional Hamiltonian systems. These findings suggest that the presence of 
chains of islands of stability in phase space and the associated diffusion and stickiness phenomena around 
these islands are responsible for QSS approximated by g-Gaussian pdfs for long times. 

Nevertheless, as these authors emphasize, it is important to recall that g-Gaussians are not the only 
possible choice for describing QSS of the type we have discussed. Indeed, in certain cases, it has been found 
that metastable states start as g-Gaussians and pass through intermediate stages where they are better 
approximated by other types of distributions befor e finally tending to Gaussians, after long enough times 
[An tonopoulos et~al. . 20 id : Ruiz-Lopez et al . 2010l ]. 

We believe, therefore, that the realm of conservative systems, represented by Hamiltonian flows or sym- 
plectic maps, is well suited for discovering complex metastable phenomena occurring in "weakly chaotic" 
domains and relating to more global properties of the dynamics. Even if the system is multi-dimensional, 
the dynamics near the boundaries of resonance islands frequently yield long-lived states with Non-Extensive 
statistics. Thus, we suggest that it may be possible to apply t hese ideas to some very slow diffusive phe- 
nomena recently o bserved in 1-dimensional disordered chains Flach et al. . 2009 : Johansson et all . 20091 : 
Skokos et ai] . l2009l ]. where it is not yet known if, in the limit of i — )■ oo, diffusion will extend to particles 



further and further away or will b e hindered by the boundary of some high-dimensional KAM torus, which 
is believed to exist [Aubrvl . [2010l ]. 



6. Appendix 



In j Christodoulidi et ad . 1201 11 ]. it is shown using the Poincare - Lindstedt perturbation method, that one can 
construct various quasi-periodic orbits for the FPU system, under the assumption that at zeroth order we 
excite s arbitrarily selected modes, with 1 < s < A^. Actually, it is possible to prove a proposition predicting 
which modes will be subsequently excited in the next orders of the perturbation theory. According to this 
proposition, if one excites in the zeroth or der of the FPU-a model (llOp the mode ko, so that s = 1 
(corresponding to the q-breather solution of [f^.c^ et all Bj) , the subsequent ki modes, where i denotes 
the order of the Poincare - Lindstedt series, will be 



ki 



{i + l)fco + N 
2{N + 1) 



(AT + 1) _ (i + i)k, 



(16) 
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where [•] denotes the integer part of the argument and the | • | the absolute value of the argument. Thus, 
for = 31 particles, one finds the sequence reported in Sec. 14.21 for the excitable /cq = 31-st mode. 
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